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INVERSION  OF  THE  ACOUSTIC  PLANE  WAVE 
REFLECTION  RESPONSE  OF  A  LAYERED  OCEAN  80TT0M 


1.  INTRODUCTION 


Proper  interpretation  of  acoustic  signals  measured  at  sea  often  requires 
that  their  interaction  within  the  ocean  bottom  he  considered.  This  is 
particularly  true  at  the  low  frequencies  used  in  passive  sonar  applications 
since  long  wavelength  sound  waves  can  penetrate  the  surficial  sediments.  The 
effect  of  a  layered  ocean  bottom  on  acoustic  propagation  is  usually  treated  in 
terms  of  bottom  loss,  defined  as  the  ratio  of  reflected  to  incident  plane  wave 
intensities  expressed  in  decibels,  as  a  function  of  frequency  and  grazing 
angle. 

Estimates  of  bottom  structure  (e.g.,  density  and  sound  speed  profiles) 
can  be  obtained  by  comparing  the  bottom  loss  inferred  from  experiments  with 
that  calculated  for  idealized  geoacoustic  models.  Parameters  of  the 
geoacoustic  models  are  adjusted  until  calculated  values  of  bottom  loss  agree 
with  measured  values  according  to  criteria  specified  for  this  curve-fitting 
approach.  Experimental  bottom  loss  is  often  inferred  from  propagation 
measuronents  in  which  small  explosive  charges  are  used  as  sound  sources. 
Although  explosives  provide  the  required  broadband  frequency  response,  their 
point-source  distribution  can  give  rise  to  multipaths  in  the  subbottom  that 
complicate  determination  of  the  bottom  loss.  Difficulties  in  interpretation 
occur  at  those  source/receiver  offsets  for  which  subbottom  refracted  arrivals 
become  time-coincident  with  arrivals  reflected  from  the  water/sediment 
interface. 1 

The  mathematical  formulation  whereby  bottom  loss  is  determined  from  a 
specification  of  the  layered  bottom  structure  is  a  well-posed  problem  of  the 
"forward"  type.  The  aim  of  investigations  of  bottom  loss  is  a  solution  of  the 
"inverse"  problem,  whereby  the  structure  of  the  ocean  bottom  is  determined 
from  a  limited  knowledge  of  the  bottom  loss.  Much  of  the  theoretical  effort 
on  the  inverse  problem  has  been  confined  to  the  model  of  plane  waves  scattered 
by  a  one-dimensional  inhomogeneous  medium. 

Two  comprehensive  reviews  of  formal  results  on  one-dimensional  inverse 
theory  are  given  by  8urridge2  and  Newton. 3  Most  of  the  effort  has 
concentrated  on  exact  methods  of  inversion,  and  few  numerical  results  on 
simulated  or  real  data  have  appeared.  An  exception  is  found  in  the  recent 
work  of  Candel  et  al.4,5  In  their  method,  the  acoustic  field  is  first  split 
into  upgoing  and  downgoing  waves. 6  Application  of  the  forward  scattering 
approximation  leads,  in  a  straightforward  way,  to  an  analytical  expression  for 
the  reflection  coefficient  in  the  form  of  a  nonlinear  Fourier  transform  of  the 
logarithmic  derivative  of  the  admittance.  Inversion  of  the  integral  transform 
leads  to  a  direct  (noniterative)  numerical  algorithm  for  determining 
admittance  as  a  function  of  depth  from  the  impulse  response  of  the  medium. 

In  this  report,  we  present  the  results  to  date  on  our  implementation  of 
the  inversion  technique  proposed  by  Candel  et  al.4*5  for  simulated  data. 
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Other  promising  numerical  inversion  schemes  that  have  appeared  in  the  recent 
1 iterature^-12  are  not  discussed.  In  section  2  we  review  Candel  et  al.'s 
treatment  of  the  forward  problem  for  global  and  local  wave  representations  and 
develop  an  equivalent  integral  equation  form.  In  section  3  we  develop  an 
analytic  expression  for  the  reflection  coefficient  based  on  the  forward 
scattering  approximation.  The  same  result  is  obtained  in  section  4  from  WKBJ 
theory.  A  brief  development  of  the  Ricatti  equation  satisfied  by  the 
reflection  coefficient  is  presented  in  the  appendix.  In  section  5  we  indicate 
how  the  approximate  expression  for  the  reflection  coefficient  leads  to  an 
efficient  inversion  method.  After  some  discussion  on  implementation  of  the 
inversion  algorithm,  we  present  several  numerical  examples  on  simulated  data 
to  illustrate  the  method  in  section  6. 


2.  THE  FORWARD  PROBLEM 


2.1  THE  MATHEMATICAL  MODEL 


The  mathematical  model  to  be  considered  is  shown  in  figure  1.  An 
inhomogeneous  liquid  with  density,  p(z),  and  sound  speed,  c(z),  which  are 
arbitrary  functions  of  depth,  occupies  the  region  0  <  z  <  H  between  two 
homogeneous  liquids.  It  is  convenient  to  regard  the  region  0  <  z  <  H  as 
divided  into  M  layers  of  thickness,  hm,  m  =  1,  2,  ...,  M.  The  density  and 
sound  speed  within  each  layer  are  taken  to  vary  continuously  with 
discontinuities  in  p{z)  and  c(z)  introduced  at  layer  boundaries.  The 
homogeneous  regions  are  characterized  by  constant  density  and  sound  speed 
pairs  p0,  c0  for  z  <  0  and  pj,  ci  for  z  >  H.  All  regions  are  assumed 
to  be  nonabsorbing. 

An  acoustic  plane  wave,  p-j,  initially  propagating  downward  at  a  grazing 
angle  e0  to  the  z  =  0  plane,  encounters  the  inhomogeneous  medium  at  z  =  0. 
Within  0  <  z  <  H,  the  incident  wave  is  partially  reflected  and  a  reflected 
wave  Rpy.  is  returned  to  z  <  0.  For  z  >  H,  where  no  further  reflections 
occur,  only  a  transmitted  wave  Tpt  continues  to  propagate. 


2.2  GLOBAL  FIELD  EQUATIONS 


For  stratified  media,  the  wave  vectors  are  confined  to,  say,  the 
xz-plane,  so  that  all  field  quantities  are  independent  of  y.  Moreover,  all 
waves  exhibit  a  common  wavenumber  in  the  x-direction,  i.e.,  kx  =  k0  cos  Op 
fixed  by  the  angle  of  the  incident  wave.  Assuming  time-harmonic  waves  wit 
angular  frequency  u  =  2»f,  the  pressure  and  particle  velocity  can  be 
represented  by  p(z)exp(ikxx  -  iwt)  and  [u(z),  0,  w(z)]exp(ikxx  -  i<ot), 
respectively.  It  follows  that  the  basic  equations  describing  the  acoustic 
field  in  an  inhomogeneous  mediuml3  1  duce  to  the  form. 


3- 
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(1) 

(2) 

(3) 


p‘  =  iwpw, 

W  =  ikzYzp, 
u  =  kxp /  (u)o ) , 

where  the  prime  denotes  d/dz  and  k,  =  (k^  -  k^)^  is  the  longitudinal 

component  of  the  total  wavenumber  k  =  w/c  at  a  given  depth  z.  Yz  defines 
a  "longitudinal  admittance," 


^z  =  (kz /k)/(oc)  -  kz/(up) .  (4) 

For  a  plane  wave  propagating  in  a  uniform  medium  with  (real)  wavenumbers  kx 
and  kz,  Yz  is  the  ratio  w/p  for  a  downgoing  wave  and  -w/p  for  an  upgoing 
wave.  At  normal  incidence  (e0  =  90°),  kz  =  k  and  Yz  reduces  to  Y  =  l/(pc). 
Equations  (1)  and  (2)  can  be  cast  into  the  compact  matrix  form. 


We  consider  only  the  case  for  which  kz  remains  real  everywhere,  i.e., 
no  total  reflection  occurs  within  the  inhomogeneous  region.  Then  the  only 
wave  propagating  in  z  >  H  can  be  represented  by 

P ( z )  =  P  exp[ikz(z  -  H) ],  (6) 

w(z)  =  W  exp[ikz(z  -  H)].  (7) 

From  equation  (2)  it  follows  that  ikzW  =  ik?YzP,  whence  W  =  YZP.  If  we 
normalize  P  to  unity,  then  appropriate  initial  conditions  for  the  system 
in  equation  (5)  are 


P(H)  =  1, 

(8) 

w(H)  =  YZ(H). 

(9) 

With  the  above  initial  conditions,  p(z)  and  w(z)  can  be  calculated 
everywhere  within  0  <  z  <  H  by  numerically  integrating  the  system  in  equation 


(5).  The  reflection  and  transmission  coefficients  may  be  obtained  by  letting 
A  and  B  designate  the  incident  and  reflected  amplitudes  in  z  <  0.  Then  the 
field  in  this  homogeneous  region  is  represented  by 

p(z)  =  A  exp(ikzz)  +  B  exp(-ikzz),  (10) 

2  2  1/2 

where  kz  =  (k0  -  kx)  .  According  to  equation  (2), 

w(z)  =  YZ[A  exp(ikzz)  -  B  exp(-ikzz)].  (11) 

Solving  equations  (10)  and  (11)  for  A  and  B  and  evaluating  at  z  =  0  give 

A  =  [p ( 0 )  +  w(0)/Yz(0)]/2,  (12) 

B  =  [p(0)  -  w(0)/Yz(0)]/2.  (13) 
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The  reflection  and  transmission  coefficients  are  then  determined  by 

R  =  B/A  =  [p(0)  -  w(0)/Yz(0)]/[p(0)  ♦  w(0) / Yz ^ 0 ) ],  (14) 

T  =  1/A  =  2/[p(0)  +  w(0)/Yz(0)].  (15) 

2.3  LOCAL  FIELD  EQUATIONS 


The  system  in  equation  (5)  together  with  initial  conditions  (8)  and  (9) 
provides  a  formal  solution  to  the  global  evolution  of  the  field.  It  is  useful 
to  obtain  an  alternate  representation  in  terms  of  local  waves,  i.e.,  at  each 
depth  z  to  decompose  the  total  field  into  upgoing  and  downgoing  components. 
Although  such  a  decomposition  is  not  unique  for  inhomogeneous  media, 14,15 
conditions  for  which  reflected  wave  amplitudes  are  small  relative  to  the 
incident  wave  amplitude  suggest  splitting  the  field  into  local  waves  as  if  the 
medium  were  locally  homogeneous. 

At  each  depth,  we  introduce  a  local  upgoing  wave,  U,  and  a  local 
downgoing  wave,  D,  defined  by 


D  =  (p  +  w/Yz)/2,  (16) 

U  =  (P  -  w/Yz)/2.  (17) 

While  the  above  decomposition  is  somewhat  arbitrary  for  inhomogeneous  media, 
it  is  a  classical  one. 6, 13, 16, 17  por  homogeneous  media,  the  splitting 
defined  by  equations  (16)  and  (17)  provides  the  desired  identification  of  two 
elementary  waves,  as  indicated  by  substitution  into  equations  (10)  and  (11). 


Equations  (16)  and  (17)  are  readily  inverted  to  give 


p  =  0  +  u. 

(18) 

w  =  Yz(0  -  U). 

(19) 

Substitution  of  equations  (18)  and  (19)  into  the  coupled  system  (5)  determines 
a  new  coupled  system  for  the  local  fields  U  and  0  in  the  form 


where  we  have  set 

9  -  Y'/(2Yz)  «  (fcnYz)'/2.  (21) 


When  the  medium  is  homogeneous,  the  system  in  equation  (20)  decouples  and 
U  and  0  take  the  usual  form  of  upgoing  and  downgoing  waves.  When  the  medium 
is  inhomogeneous,  coupling  between  U  and  D  arises  via  g,  the  relative  variation 
of  the  longitudinal  admittance. 
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Integration  of  the  system  in  equation  (20)  can  be  started  at  z  =  H. 

Below  that  depth,  the  only  wave  is  downgoing.  Choosing  the  downgoing 
amplitude  as  unity  determines  the  initial  conditions, 

0(H)  =  1,  (22) 

U(H)  =  0,  (23) 

for  the  system  in  equation  (20),  which  accord  with  conditions  (8)  and  (9)  for 
the  system  in  equation  (5)  obtained  earlier.  Since  U(z)  and  D(z)  are 
determined  at  every  depth  z  by  numerical  integration,  it  is  reasonable  to 
define  local  reflection  and  transmission  coefficients  as 

R(z)  =  U(z)/D(z),  (24) 

T(z)  =  1/D(z) .  (25) 

The  global  field  quantities  p  and  w  are  continuous  across  discontinuities 
in  the  properties  of  the  medium.  The  local  waves  U  and  D,  however,  are 
discontinuous  there.  Accordingly,  some  care  is  required  when  discontinuities 
in  density  and/or  sound  speed  are  encountered  during  the  integration  of  the 
system  in  equation  (20).  Suppose  an  admittance  jump  of  magnitude  Y+  -  Y_ 
occurs  at  depth  zQ,  where  Y±  =  Y^(zo  ±  0).  Then  the  continuity  of  p  and~w 

together  with  equations  (18)  and  (19)  determines  the  local  wave  U_  and 
D  at  z„  -  0  in  terms  of  the  local  waves  U.  and  D,  at  z  +  0  according  to 

P_  =  0_  +  U_  =  D+  +  U+  -  p+,  (26) 

w_  =  Y_(D_  -  UJ  =  Y+(D+  -  U+)  =  w+,  (27) 

from  which  the  required  jump  conditions  on  the  local  waves  are  found  to  be 

0_  =  0+(l  +  Y+/Y_)/2  +  U+(l  -  Y+/Yj/2,  (28) 

U_  =  D+(l  -  Y+/Y_)/2  +  U+(l  ♦  Y+/Y_)/2.  (29) 


2.4  INTEGRAL  EQUATION  REPRESENTATION 


The  system  in  equation  (20)  can  be  rearranged  into  the  alternate  form. 
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The  differential  system  in  equation  (30)  together  with  initial  conditions  in 
equations  (22)  and  (23)  can  be  converted  into  an  equivalent  integral  form, 
e.g.,  by  the  variation  of  parameters  method,  to  obtain 

D(z)  =  exp[f+(z)]  (1  +  /  g(s)  U(s)  exp[-f+(s)])  ds,  (32) 

H 

U(z)  =  exp[f_(z)]  /  g(s)  D(s)  exp[-f_(s)]  ds,  (33) 

H 

where 


f±(z) 


z 

/g±(s)  ds. 


H 


(34) 


To  this  point  the  theory  is  exact.  Approximations  to  U  and  D  are  now 
developed,  leading  to  a  useful  analytic  result  for  R,  which  forms  the  basis  of 
a  noniterative  method  of  profile  inversion.  It  is  worthwhile  noting  that, 
whereas  the  local  waves  U  and  D  satisfy  the  differential  system  in  equation 
(20),  the  reflection  coefficient  R  =  U/D  satisfies  a  nonlinear  equation  of  the 
Ricatti  type.  This  point  is  developed  briefly  in  the  appendix. 
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3.  FORWARD  SCATTERING  APPROXIMATION 


While  the  differential  system  in  equation  (20)  is  readily  solved 
numerically,  an  interesting  approximation  method  may  he  developed  that  leads, 
in  a  straightforward  way,  to  an  analytic  representation  for  the  reflection 
coefficient.  It  is  convenient  to  derive  this  result  hy  solving  equations  (32) 
and  (33)  by  successive  approximations.  In  this  scheme,  higher  order  iterates 
u(i+i)  and  u(i+D  are  obtained  from  previous  ones,  u(i)  and  u(i),  by 
substituting  the  latter  values  into  the  right-hand  side  of  equations  (32)  and 
(33).  Candel  et  al.4  initiate  the  procedure  by  using  a  forward  scattering 
approximation.  This  consists  of,  first,  neglecting  the  upgoing  wave  d  front 
equation  (32),  solving  for  L),  and  then  using  the  D  obtained  in  this  step  in 
equation  (33)  and  solving  for  U.  Although  it  is  clearly  possible  to  continue 
iterating  in  this  fashion,  for  |ll  j  <<  |  0|  it  is  useful  to  examine  the 
approximate  solution  obtained  at  this  stage.  In  the  approach  just  described, 
the  differential  system  in  equation  (20)  is  replaced  with  a  system  for  which 
the  upper  off-diagonal  element  -1  on  the  right-hand  side  is  replaced  by  zero. 

From  equation  (32)  with  tl(0)  =  0,  the  solution  for  i)(0)  is  found  to  be 

z 

[)(°)(z)  =  exp[f+(z) ]  =  [Yz(H)/Yz]1/2  exp[i  /  kz(s)  ds]. 

H 

Substitution  of  equation  (35)  into  the  right-hand  side  of  equation 
to  the  solution  for  llU)  in  the  form 

z  s 

l/^(z)  =  exp[f_(z)l  /  g(s)  exp[2i/kz(t)  dt]  ds.  (36) 

H  H 

A  useful  approximation  for  the  local  reflection  coefficient  can  now  be 
obtained  from  equations  (35)  and  (36)  by  forming  the  ratio  u( 1) ( z) /d(°) ( z) . 

The  result  evaluated  at  z  =  0  is 

H  s 

R^(0)  =  -/g(s)  exp[2i/kz(t)  dt]  ds,  (37) 

0  0 

where  the  limits  of  integration  have  been  altered  according  to 

s  z  s 
H  H  z 

Althouqh  equation  (37)  is  an  approximation,  it  has  been  shown 
numerically**  to  provide  accurate  results  in  many  situations  of  practical 
interest.  Moreover,  this  simple  result  obtained  from  the  forward  scattering 
approximation  forms  the  basis  of  a  method  to  recover  the  admittance  versus 
depth  of  an  inhomogeneous  medium  from  the  impulse  response.  We  will  develoD 
too  formalism  of  this  method  in  a  subsequent  section. 
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(35) 
(33)  leads 
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It  is  shown  in  the  appendix  that  equation  (37)  follows  from  the  Ricatti 
equation  formulation  for  R  =  U/D.  Before  proceeding  to  the  development  of  the 
inversion  approach  from  equation  (37),  it  is  worthwhile  to  examine  this  result 
within  the  context  of  the  WKBJ  approximation. 
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4.  WKBJ  APPROXIMATION 


The  approximation  given  by  equation  (37)  has  been  derived 
elsewhere. 8, 13, 14  it  can  be  developed  from  a  WKBJ  representation  for  the 
solution  in  the  region  0  <  z  <  H.  The  appropriate  WKBJ  form  for  the  acoustic 
problem  is 


p(z)  =  A(z)W+(z)  +  B(z)W_(z) ,  (38) 

where  A(z)  and  B ( z )  are  functions  to  be  determined  and 

W±(z)  =  exp[±  i/kz(s)  ds]/[Yz(z)]1/2.  (39) 

H 

(W±  is  not  to  be  confused  with  the  amplitude  of  the  particle  velocity  used 
in  section  2.2.)  We  recall  that  since  we  are  only  considering  real  propagating 
waves,  the  denominator  of  equation  (39)  is  well  defined  throughout  0  <  z  <  H. 

A  coupled  system  of  differential  equations  for  A  and  B  is  readily  obtained 

by  direct  manipulation.*6  With  the  identification  D  =  AW+  and  U  =  BW_  ,  there 
follows  from  equations  (30),  (31),  (38),  and  (39)  the  sequence 

O'  =  (-g  +  ik2)D  +  gU 

=  (W;/W+)D  +  gU 

=  aw;  +  gBW_  =  A'W+  +  AW|, 


which  determines  the  result 


A*  =  gBYzW^.  (40) 

The  result  for  B1  is  obtained  in  a  similar  way.  The  final  coupled  system 
for  A  and  B  can  be  expressed  in  the  compact  matrix  form. 


From  the  initial  conditions  for  U  and  D,  it  follows  that  appropriate  initial 
conditions  for  the  system  in  equation  (41)  are  given  by 


A(H)  =  [Yz(H)]1/2. 

(42) 

• 

o 

II 

00 

(43) 
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The  above  initial  conditions  are  readily  incorporated  into  the  integral 
equation  form  equivalent  to  the  differential  system  in  equation  (41),  i.e., 

Mz)  =  [YZ(H)]1/2  +  ;g(s)  B(s)  exp[-2i  /  kz(t)  dt]  ds,  (44) 

H  H 

z  s 

B(z)  =  /g(s)  A(s)  exp[+2i/k  (t)  dt]  ds,  (45) 

H  H  z 

The  usual  assumption  invoked  in  attempting  a  solution  of  the  WKBJ  form 
given  by  equations  (38)  and  (39)  is  that  variations  in  the  medium  are  slow 
compared  with  the  wavelength  of  the  waves.  Under  these  conditions,  the 
relative  variation  of  longitudinal  admittance  |g  |  «  1  and  the  system  in 
equations  (44)  and  (45)  can  be  solved  by  successive  approximation.  Setting 
g  =  0  initially  leads  to  zeroth  order  estimates  A(0)  =  [yz(h)]1/2  and 
b(0)  *  0.  Substituting  these  results  into  the  right-hand  sides  of  equations 
(44)  and  (45),  we  find  AU)  =  a(0)  and 

B^(z)  =  [Yz(H)]1/2  /  g(s)  exp[2i/kz(t)  dt]  ds.  (46) 

H  H 

The  approximate  first  order  reflection  coefficient  given  by  R ( 1 )  = 
(B^W_)/(A^W+)  when  evaluated  at  the  surface  z  =  0  takes  the  form 


R  1  (0)  =  -fg(s)  exp[2i  f  k  (t)  dt]  ds. 


which  agrees  with  the  result  obtained  previously  via  the  forward  scattering 
approximation. 

It  is  worthwhile  remarking  that  the  system  in  equation  (41)  can  be  solved 
numerically.  McKisic  and  HammlS  suggest  one  method.  For  the  problems  to  be 
considered  in  the  present  work,  with  kz  everywhere  real,  a  simple 

1/2 

formulation  is  obtained  by  setting  C  =  Y  W+,  from  which  we  obtain  the  system 


'0  g/c2 

gc2  o 


o  I  B 


subject  to  initial  conditions 


A(H)  =  [Yz(H)]1/2, 
B(H)  =  0, 


C(H)  =  1. 
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Jump  conditions  on  A  and  B  must  be  applied  whenever  discontinuities 
in  Yz  are  encountered.  From  equations  (26)  and  (27)  and  the  relationships 

0  =  AC/Y*/2  and  U  =  B/(Cy]/2),  the  required  connections  are  readily  deduced  to 

be 

A_  =  [Y_/Y+]1/2  [A+( 1  +  Y+/Y_)/2  +  C-2  B+(l  -  Y+/Y_)/2],  (52) 

B_  =  [WY+]1/2  [C2  A+( 1  -  Y+/Yj/2  +  B+(l  +  Y+/Y_)/2],  (53) 


where  C  =  C  =  C,. 
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5.  THE  INVERSE  SOLUTION 


5.1  PHYSICAL  INTERPRETATION 


Before  proceeding  with  the  development  of  the  inversion  scheme  based  on 
equation  (37),  it  is  worthwhile  commenting  on  the  physical  interpretation  of 
this  approximate  result.  Since  all  subsequent  work  will  focus  on  the 
reflection  coefficient  at  z  =  0,  we  drop  the  depth  dependence  from  the 
argument  of  R  and  replace  it  with  the  frequency  dependence  since  it  is 
apparent  that  R  depends  on  f  via  the  parameter  kz. 

The  reflection  coefficient  in  equation  (37)  appears  as  a  sum  of 
contributions  from  successive  elementary  layers.  About  the  layer  at  depth 
z  =  z0,  we  obtain  a  partial  reflection  coefficient, 

2o 

dR(f)  =  -  g(z  )  exp[2i  / k  (s)  ds]  dz  , 
o  o  2  o 

with  an  amplitude  determined  by  the  relative  variation  of  longitudinal 
admittance  at  that  depth  and  a  phase. 


o 

d(f)  =  2s  k  (s)  ds, 

0  z 

associated  with  the  time  taken  for  the  incident  wave  to  travel  to  depth  z  =  z0 
and  to  return  to  depth  z  =  0  as  a  partially  reflected  wave.  This  is  seen  more 
clearly  by  writing  tf(f)  =  2irf-r,  where 


t  =  2^  ds/c  (s)  (54) 

0  z 


is  the  total  propagation  time  from  z  =  0  to  z  =  z0  and  back  for  a  wave 
having  longitudinal  phase  speed  c2  =  a>/kz. 


5.2  THE  INVERSION  FORMULAS 


We  observe  that  equation  (37)  determines  the  reflection  coefficient  as  a 
nonlinear  Fourier  transform.  Inverting  this  Fourier  transform  recovers  the 
relative  variation  of  the  longitudinal  admittance  and  integration  determines 
the  admittance  itself.  We  remark  here  that  only  Yz  can  be  recovered  from  a 
knowledge  of  R(f)  for  a  given  angle  of  incidence.  Independent  determination 
of  both  the  density,  p,  and  the  sound  speed,  c,  requires  more  information.  We 
will  return  to  this  point  later. 
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To  cast  equation  (37)  in  the  form  of  a  standard  Fourier  transform,  we 
introduce  a  new  coordinate  t,  defined  by 

C  =  (2/ko) /kz(s)  ds,  (55) 

where  k0  designates  the  wavenumber  in  region  z  <  0.  With  the  above  change 
of  variable,  we  obtain  the  result 


C(H) 

R(f)  =  -  /  g(c)  exp(2irif?/c  )  dc,  (56) 

0  0 

where  g(c)  =  Y^(£ ) /[2Yz(? ) ]  and  the  prime  is  now  used  to  denote  d/dc.  The 

reflection  coefficient  now  appears  as  a  standard  Fourier  transform  of  the 
relative  variation  of  longitudinal  admittance  with  respect  to  the  new  variable 

t. 


To  invert  equation  (56)  we  first  compute  the  time  response  of  the  medium, 
i.e.,  the  reflection  coefficient  as  a  function  of  time  given  by 

oo 

r(t)  =  /  R(f )  exp(-2irift)  df.  (57) 

-oo 

Transforming  both  sides  of  equation  (56)  leads  to  the  result 

£  ( H)  <x> 

r(t)  =  -  I  I  gk)  exp[2irif(c  /C  -  t]  df  d£.  (58) 

0  -«  0 

The  integration  with  respect  to  f  may  be  performed  at  once  giving 

oo 

/  exp[2nifU/c0  -  t)]  df  =  c0«U  -  c0t),  (59) 

—  00 

whence  the  sifting  property  of  the  delta  function  determines  the  time  response 
in  the  form 

r(t)  =  -  c0g(c0t)  =  -  c0[(dYz/dC)/(2Yz)]  |.  (60) 

C  =  c0t 

The  physical  interpretation  of  equation  (60)  indicates  that  the  time 
response  at  instant  t  is  proportional  to  the  relative  variation  of  admittance 
at  the  layer  of  coordinate  C  =  c0t  corresponding  to  an  actual  location  given 
by 

cot  =  (2/k0) /  kz(s)  ds  (61) 


or 


t  =  2 /  ds/c  (s) . 


(62) 
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The  "active"  layer  at  instant  t  is  located  at  a  distance  z  corresponding  to 
the  travel  time  from  z  =  0  to  this  layer  and  back. 

Equation  (60)  may  be  rearranged  in  the  form 

Vznz  =  -2r(t)/c0,  (62 

S=cot 

and  equation  (55)  may  be  replaced  by 

z'  =  k0/(2kz)  =  l/(2pc0Yz),  (64 

where  use  was  made  of  equation  (4).  Since  these  two  relations  (with 
appropriate  initial  conditions)  allow  the  determination  of  Yz  versus  z,  the 
direct  inversion  scheme  is  formally  complete. 

The  inversion  algorithm  of  Candel  et  al.5  is  formed  from  equations 
(4),  (63),  and  (64).  Three  cases  have  to  be  considered  separately. 


5.2.1  Case  A;  p(z)  Known,  c(z)  Unknown 


When  density  p(z)  is  known,  the  sound  speed  c(z)  can  be  determined  by 
integrating  equations  (63)  and  (64)  directly,  i.e., 

Y;  =  -(2/c0)rYz,  (6! 

z'  =  l/(2pC0Yz) .  (6( 

At  each  step  of  the  numerical  integration,  the  sound  speed  can  be  recovered 
from  equation  (4),  which  can  be  written  in  the  convenient  form 

n  =  02coYz  +  cos2e0]1/2,  (67 

where  n(z)  =  c0/c(z)  =  k(z)/k0  is  the  local  index  of  refraction  and 
cos  e0  =>  kx/k0  determines  the  grazing  angle  at  z  =  0.  For  a  receiver 
located  at  S0,  corresponding  to  a  point  at  z  =  z0  in  z  <  0,  the 
integration  of  equations  (651  and  (66)  may  be  started  with  the  initial 
conditions 

Yz(S0)  =  sin  o0/(poco).  (6? 

zo  =  40/(2  sin  e0 ) .  (65 


5.2.2  Case  B:  c(z)  Known,  p(z)  Unknown 


For  the  case  when  the  sound  speed,  c(z),  is  known  and  the  density,  p(z). 
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is  to  be  determined,  equations  (4)  and  (64)  have  to  be  modified  slightly,  but 


equation  (63)  remains  unchanged. 

We  find 

Y*  = 
z 

-(2/c0)rY2, 

(70) 

z 1  = 

[n2  -  cos2 

9 0]"1/2/2’ 

(71) 

p  = 

[n2  -  cos2 

9o]1/2/(Vz)- 

(72) 

The  system  in  equations  (70) 

and  (71)  is 

directly  integrable 

since  n(z) 

is  known  everywhere.  At  each  step  in  the  integration,  the  density  is  obtained 
simply  from  equation  (72).  The  initial  conditions  for  equations  (70)  and  (71) 
are  the  same  as  in  case  A,  subsection  5.2.1. 


5.2.3  Case  C:  p(z)  and  c(z)  Unknown 


As  indicated  previously,  unique  determination  of  both  density,  p(z),  and 
sound  speed,  c(z),  requires  additional  information.  For  the  application  of 
Candel  et  al.'s^  inversion  method,  reflection  responses  for  at  least  two 
distinct  grazing  angles  must  be  provided. 

Let  the  subscripts  1  and  2  denote  the  quantities  that  correspond  to  the 
two  distinct  grazing  angles  ©i  and  ©2,  with  »2  >  el*  The  relevant 
differential  equation  system  for  grazing  angle  ©^  is  given  by 

(Yz){  -  -(2/c0)r1(Yz)1,  (73) 


z*  =  l/[2Pco(Yz)1], 


(74) 


(Yz)x  =  [n2  -  cos2  e1]1/2/(pco). 


(75) 


where  the  prime  now  denotes  d/d^i .  A  similar  set  of  equations  corresponds 
to  the  other  grazing  angle  ©2,  but  it  is  apparent  that  the  integration 
variables  ci  and  t2  are  n°t  identical.  Use  of  the  chain  rule  for 
differentiation,  however,  relates  to  ?2  so  that  the  system  for  grazing 
angle  ©2  can  be  specified  in  terms  of  ci  in  the  form 

r’2  =  (Yz)2/(Yz)1’  (76) 

( yz}2  =  -(2/cQ)r2(Yz)|/( Yz)1,  (77) 


(Yjo  =  [n2  -  cos2  ©o]^2/(pcJ. 


(78) 
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Now  is  determined  as  a  depenuent  variable  that  allows  (Yz)i  and 

( Yz ) 2  to  be  determined  at  the  same  depth  location  z.  The  density  and 
sound  speed  are  then  readily  deduced  by  combining  equations  (75)  and  (78): 

p2c2  =  [cos2  -  cos2  e2]/[(Yz)2  -  (Yz)2],  (79) 

n2  =  [(Yz)2  cos2  e1  -  (Yz)2  cos2  &2]/[(Yz)2  -  (Yz)2].  (80) 

It  is  evident  that  the  integration  of  equation  (74)  can  proceed  only  in 
conjunction  with  the  determination  of  p  via  equation  (79).  Appropriate 


initial  conditions  for  this  case  become 

Yz(?0)l  =  sin  ei/(p0c0).  (81) 

zo  =  W(2  sin  ®l)»  (82) 

^2(^o)  =  ?oYz(co)2/Yz(^o)l»  (83) 

yz(^o)2  =  sin  ©2/(poCq)*  (84) 


In  principle,  it  is  no  more  difficult  to  determine  both  the  density, 
p(z),  and  the  sound  speed,  c(z),  profiles  than  one  of  them.  However, 
reflection  responses  for  at  least  two  probing  directions  are  required,  and 
four  instead  of  two  differential  equations  must  be  integrated. 


5.3  REMARKS 


Before  proceeding  with  numerical  aspects  of  both  the  forward  and  inverse 
problems,  some  comments  on  the  analysis  so  far  are  in  order. 

The  steps  leading  to  equation  (6)  required  the  relative  variation  of 
admittance  to  be  independent  of  frequency.  While  this  is  most  certainly  a 
valid  assumption  for  the  density  and  sound  speed,  usual  treatments  regarding 
absorption  presume  at  least  a  linear  dependence  on  f.  Therefore,  absorption 
was  taken  to  be  zero  in  the  above  analysis.  It  may  be  possible  to  incorporate 
some  frequency  dependence  into  Yz  and  still  manipulate  the  integrands  in  a 
straightforward  way,  but  that  task  will  not  be  undertaken  in  this  report.* 
Another  reason  for  neglecting  the  effects  of  absorption  lies  in  numerical 
difficulties  associated  with  integrating  differential  equations  containing  a 
damping  term. 

The  restriction  to  real  propagating  waves  everywhere  within  the 
inhomogeneous  region  0  <  z  <  H  is  partly  due  to  numerical  difficulties 
associated  with  turning  points  within  the  domain  of  integration.  Ad  hoc 
devices  are  sometimes  introduced  to  avoid  numerical  underflow/overflow 

*  In  fact,  if  absorption  in  the  subbottom  is  proportional  to  frequency,  then 
Yz  remains  independent  of  f. 


problems  in  these  situations. 19  On  the  other  hand,  waves  are  exponentially 
attenuated  below  a  turning  point  and  quickly  lose  any  practical  capability  to 
return  information  from  regions  at  greater  depths. 

In  the  WKBJ  approach  to  the  approximate  analytical  representation  for  the 
reflection  coefficient,  the  criterion  |g  |  «  1  was  used.  This  condition  seems 
to  preclude  regions  containing  discontinuities  of  the  first  kind  from 
consideration.  However,  the  forward  scattering  approach  of  Candel  et  al.4 
only  assumes  |U|  «  |D|,  which  can  hold  even  in  the  presence  of  discontinuous 
jumps  in  material  properties.  Clearly  the  latter  condition  encompasses  the 
former  one.  It  is  interesting  that  both  criteria  lead  to  the  same  approximate 
formula  for  the  reflection  coefficient. 

The  inversion  scheme  suggested  by  equations  (63)  and  (64),  which  only 
requires  integration  of  a  pair  of  first  order  differential  equations,  together 
with  equation  (4),  is  direct  since  its  implementation  does  not  involve  an 
iteration  procedure.  On  the  other  hand,  it  would  appear  that  the  scheme  only 
accounts  for  first  order  reflections  from  the  inhomogeneous  medium.  This  is 
partly  because  the  mapping  between  the  depth  coordinate,  s,  and  the  time,  t, 
is  monotonic;  i.e.,  later  times  correspond  to  deeper  depths.  For  conditions 
in  the  subbottom,  which  support  multiple  scattering,  the  correspondence  is  not 
monotonic.  Because  of  this  limitation,  configurations  of  material  properties 
giving  rise  to  multiple  scattering  should  be  monitored  during  the  inversion 
process  and  their  effect  taken  into  account  in  the  interpretation  of  results. 
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6.  NUMERICAL  RESULTS 


In  this  section,  we  present  some  numerical  results  based  on  a  computer 
implementation  of  Candel  et  a  1 .  ‘ s5  inversion  procedure  for  simulated  data. 
Their  numerical  examples  were  restricted  to  waves  propagating  at  normal 
incidence  (e0  =  90°)  in  regions  of  constant  density  (o(z)  =  p<0  only,  and 
recovery  of  n(z)  required  only  a  single  impulse  response,  r(t).  Here  we 
extend  the  numerical  results  to  incorporate  two  impulse  responses 
corresponding  to  both  normal  and  oblique  grazing  angles  so  that  simultaneous 
recovery  of  both  n(z)  and  p(z)  is  possible.  For  most  applications  to  acoustic 
probing  of  the  ocean  bottom,  both  n(z)  and  p(z)  are  usually  unknown. 

For  each  model  to  be  considered,  the  computer  simulation  is  carried  out 
in  the  following  way.  At  a  given  grazing  angle,  the  bandlimited  frequency 
response  R ( f )  is  determined  for  a  finite  set  of  discrete  frequencies 
fk  =  kAf,  k  =  0,  1,  ...,  K-l.  The  reflected  time  series  r(tn),  tn  =  nat, 
n  =  0,  1,  ...,  N-l  is  determined  via  a  discrete  inverse  Fourier  transform. 

The  discrete  transform  is  effected  using  a  fast  Fourier  transform  (FFT) 
algorithm.  N-K  zeros  are  appended  to  each  R ( f ^ )  response  to  enhance  the 
resolution  of  r(t^)  due  to  the  FFT  constraint,  Ataf  =  1/N.  With  rj  and 
r2  computed  for  distinct  grazing  angles  ej  and  03,  with  >  el* 
recovery  of  n(z)  and  p(z)  proceeds  according  to  case  C,  subsection  5.2.3. 
Reconstructed  n(z)  and  p(z)  are  finally  compared  with  the  refractive  index  and 
density  profiles  initially  used  to  generate  the  synthetic  responses. 

For  the  forward  problem,  any  of  the  four  representations  of  the  acoustic 
field  (system  in  equation  (5)  for  p  and  w;  system  in  equation  (20)  for  D  and 
U;  system  in  equation  (48)  for  A,  B,  and  C;  and  equation  (A-l)  for  R)  together 
with  appropriate  initial  and  jump  conditions  can  be  used  to  determine  the 
reflection  coefficient.  Experience  with  each  system  indicates  that  the 
representation  in  terms  of  the  local  waves  D  and  U  requires  the  least 
computational  effort.  For  both  the  forward  problem  and  its  inverse  solution, 
numerical  integration  of  the  relevant  system  of  first  order  differential 
equations  was  carried  out  using  a  computer  code  developed  by  Shampine  and 
Gordon. 20  This  cocje  is  well  documented  and  features  automatic  local  error 
control.  Moreover,  it  has  been  used  successfully  in  other  studies  of  plane 
wave  reflection  coefficients  for  layered  geoacoustic  models. 21  To 
accommodate  complex  functions  such  as  D  and  U  appearing  in  the  system  in 
equation  (20),  it  is  first  necessary  to  separate  the  equations  into  real  and 
imaginary  parts.  In  addition,  for  numerical  work,  it  is  recommended  that  all 
equations  be  converted  to  a  dimensionless  form.  This  is  easily  accomplished 
by  introducing  the  scaled  quantities  z*  =  z/H,  f*  =  fH/c0,  and  y£  =  p0c0Y7, 
where  the  scaling  factors  are  the  thickness  of  the  inhomogeneous  region,  R, 
and  the  density,  p0,  and  sound  speed,  c0,  of  the  region  z  <  0.  For  the 

inverse  problem,  we  also  have  ?*  =  ?/H  =  c0t/H  =  t*.  As  a  result  of 

this  scaling,  it  was  possible  to  maintain  numerical  accuracy  using  single 
precision  arithmetic. 

For  the  results  presented  in  this  section,  the  sound  speed  and  density 
profiles  within  the  inhomogeneous  region  0  <  z  <  H  are  modeled  in  the 

following  way.  For  each  model,  the  region  0  <  z  <  H  is  divided  into  M  layers 
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of  thickness,  hm,  m  =  1,  2,  M.  Within  the  mth  layer,  the  density  is 
assumed  to  be  constant  while  the  variation  of  c(z)  is  modeled  according  to 
l/c(z)  =  amz  +  bm.  The  layer  coefficients  am  and  bm  are  determined 
from  the  respective  sound  speeds  at  the  top  and  bottom  of  the  mth  layer.  With 
this  prescription,  the  admittance  varies  linearly  within  each  layer. 
Discontinuities  in  c(z)  or  p(z)  can  be  introduced  at  layer  interfaces. 

Three  models  are  considered  in  the  numerical  examples.  Model  1  comprises 
a  linear  refractive  index  and  constant  density  layer.  Model  2  provides  an 
example  of  a  sudden  discontinuity  in  both  refractive  index  and  density.  Both 
of  these  models  were  discussed  by  Candel  et  al.4,5  Model  3  represents  a 
geoacoustic  model  of  Hatteras  Abyssal  Plain  based  on  traditional  seismic 
interpretation  of  time  waveforms. 22  The  sound  speed  and  density  profiles 
for  the  three  models  are  summarized  in  figure  2. 

For  each  model,  the  complex  frequency  response  R  =  Re[R]  +  i  Im[R]  was 
computed  at  256  discrete  frequencies  f|<  =okfif,  k  =  0,  1,  ...,  255,  and 
Af  =  0.5  Hz  for  grazing  angles  60°  and  90°.  The  complex  sequence  R(f)  was 
extended  to  N  =  1024  points  by  appending  512  zero  values.  Since  the  time 
response  r(t)  is  a  real  sequence,  R(f)  satisfies  the  symmetry  conditions 
Re[R(-f)]  =  Re[R(f)]  and  Im[R(-f)]  =  -  Im[R{f ) ] . 5 , 23  Complex  arithmetic  was 
avoided  by  using  an  FFT  algorithm  specially  designed  to  treat  discrete 
transforms  of  real  sequences  and  their  inverses. 23,24  The  inverse  FFT 
produced  1024  estimates  of  the  time  response  at  tn  =  nAt,  n  =  0,  1,  ..., 

1023,  where  At  =  l/(Naf)  =  1/512  seconds.  Each  time  response  was  then 
convolved  with  a  low-pass  digital  filter  designed  using  the  window 
method. 25  Filters  using  both  rectangular  and  Kaiser  windows  were  employed 
with  a  cutoff  frequency  of  64  Hz.  Finally,  the  filtered  time  responses  were 
multiplied  by  Af  to  approximate  the  analytical  Fourier  transform  results. 5»23 


6.1  MODEL  1 


Figures  3  and  4  depict  the  frequency  responses  of  the  reflection 
coefficient  for  model  1.  The  amplitude  and  phase  of  R(f)  are  shown  for  a 
grazing  angle  of  60°  in  figure  3  and  90’  in  figure  4.  For  both  grazing  angles 
|R|  decreases  globally  with  increasing  frequency  and  exhibits  weak  oscillations. 
The  phase  increases  monotonical ly  but  it  is  nonlinear.  It  is  easy  to  see  that 
| R J  is  greater  for  oblique  than  for  normal  incidence.  At  zero  frequency,  the 
values  of  |  R  |  correspond  to  the  values  obtained  for  a  sudden  discontinuity  in 
the  sound  speed  of  magnitude  c0/ci.  This  example  is  discussed  in  greater 
detail  by  Candel  et  al.4,5  for  the  case  of  normal  incidence. 

The  time  responses  for  model  1  are  shown  in  figure  5.  The  upper  trace 
corresponds  to  e0  =  60“  and  the  lower  to  eg  =  90".  Each  response  has  been 
normalized  by  its  peak  value  of  r.  The  initial  time  offset  corresponds  to  the 
two-way  travel  time  for  a  receiver  located  at  z  =  -100  m.  The  time  delay 
t0  =  (2z/c0l  sin  On  is  implemented  by  applying  the  complex  modulation 
exp(-2nifT0;  to  R(f)  before  taking  the  inverse  Fourier  transform.  The 
smaller  time  delay  for  oblique  incidence  follows  from  equation  (54)  since  the 
wave  propagates  at  the  longitudinal  phase  speed  cz.  The  initial  rise  in  each 
trace  is  due  to  the  discontinuity  in  the  sound  speed  gradient  at  z  -  0. 
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Figure  5.  Bandlimited  (0-64  Hz)  Time  Responses  of  the  Reflection 
Coefficients  for  Model  1  at  Grazing  Angles  (a)  ©0  =  60*  and  (b)  e0  =  90', 

Rectangular  Window  Filter 
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Continuous  reflections  are  observed  as  the  wave  traverses  the  inhomogeneous 
layer.  The  duration  of  the  time  response  is  less  for  oblique  incidence,  but 
there  is  a  greater  peak  response  at  the  di scontinuity  in  the  sound  speed 
gradient  at  z  =  H.  The  slight  oscillations  evident  in  both  responses  are  a 
result  of  truncating  the  frequency  responses  with  a  rectangular  window  at  the 
cutoff  frequency  of  64  Hz.  The  discontinuity  in  R ( f )  generated  by  this 
filtering  operation  gives  rise  to  the  well-known  Gibbs  phenomenon. 23, 25 
These  oscillations  are  small  for  this  example  since  I R|  is  small  at  cutoff. 

The  inversion  algorithm  of  Candel  et  al.5  Was  used  to  determine  the 
reconstructed  values  of  p(z)  and  n(z)  in  figure  6  from  the  time  responses  in 
figure  5.  For  comparison,  the  density  and  sound  speed  profiles  of  the  model 
used  to  generate  R(f)  are  displayed  on  the  same  graphs.  It  is  evident  that 
the  inversion  scheme  gives  excellent  results  for  this  example. 


6,2  MODEL  2 


The  frequency  responses  for  the  step-discontinuity  profile  of  model  2  are 
shown  in  figures  7  and  8.  In  this  example,  both  density  and  sound  speed  have 
a  jump  discontinuity  at  z  =  50  m.  The  responses  for  e0  =  60'  and  »0  *  90“ 
are  given  in  figures  7  and  8,  respectively.  For  this  model,  the  analytical 
form  of  R  is  known,4  and  the  numerical  results  correctly  reproduce  the 
constant  amplitude  and  linear  phase  behavior.  The  reflection  amplitude  is 
again  greater  for  oblique  incidence.  The  normalized  time  responses  obtained 
using  a  rectangular  window  filter  are  shown  in  figure  9  for  a  receiver  located 
at  z  =  -100  m.  As  before,  the  peak  response  is  larger  and  arrives  earlier  for 
oblique  incidence.  For  this  example,  the  discontinuity  introduced  by 
truncating  R  at  64  Hz  produces  significant  Gibbs  oscillations.  Instead  of  a 
"pulse-like"  arrival,  a  sine-response^  is  observed. 

Figure  10  shows  the  reconstruction  of  o(z)  and  n(z)  based  on  the  time 
responses  in  figure  9.  Although  the  reconstructed  values  reproduce  the  glooal 
behavior  of  the  model  values,  the  effect  of  the  large  Gibbs  oscillations  is 
apparent.  Since  these  oscillations  result  from  the  discontinuity  introduced 
by  the  rectangular  filter,  it  it  reasonable  to  design  a  filter  to  reduce  this 
effect.  Figure  11  shows  time  responses  for  model  2  after  filtering  with  a 
Kaiser  window  filter.  For  these  results,  a  31-point  filter  with  60-d8 
sidelobe  suppression  was  used. 25  it  is  evident  that  most  of  the  effects  of 
the  truncation  have  been  removed.  The  reconstruction  based  on  the  Kaiser- 
windowed  time  responses  in  figure  11  are  displayed  in  figure  12.  Except  for 
the  values  near  the  interface  itself,  the  degradation  due  to  the  Gibbs 
oscillations  has  been  removed.  For  z  >  50  m,  the  reconstructed  estimate  for 
n(z)  is  slightly  in  error. 
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Figure  7.  (a)  Amplitude  and  (b)  Phase  of  the  Passband  Frequency  Response  of 

the  Reflection  Coefficient  for  Model  2,  Grazing  Angle  e0  =  60* 
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I 

6.3  MODEL  3 


The  final  example  is  based  on  a  realistic  geoacoustic  model  for  the 
Hatteras  Abyssal  Plain  region.  This  model  was  deduced  from  seismic  analysis 
of  time  waveforms  measured  at  many  source-receiver  ranges.  The  reconstruction 
below  is  based  on  impulse  responses  for  two  distinct  grazing  angles 
corresponding  to  only  two  source-receiver  offsets.  The  frequency  responses 
computed  for  this  model  are  shown  in  figure  13  for  e0  =  60*  and  in  figure  14 
for  e0  =  go*.  In  these  responses  the  modulations  are  indicative  of  multiple 
reflections.  The  nonlinear  phases  result  from  the  sound  speed  gradients 
within  the  thick  layer.  At  the  higher  frequencies  the  large  values  of  |R  [ 
suggest  the  presence  of  discontinuities  within  the  inhomogeneous  region. 

Figure  15  shows  the  time  responses  for  this  model  when  a  rectangular 
window  filter  is  used.  The  three  "pulses"  observed  at  both  grazing  angles 
correspond  to  reflections  from  discontinuities  at  z  =  0,  20.4,  and  300  m. 
Between  the  second  and  third  reflections,  continuous  returns  are  observed  from 
the  deep  layer.  The  Gibbs  oscillations  are  pronounced.  The  reconstructed 
values  of  p(z)  and  n(z)  based  on  the  time  series  in  figure  15  are  given  in 
figure  16  together  with  the  model  values.  In  spite  of  the  Gibbs  oscillations, 
the  global  behavior  of  the  reconstructions  agree  well  with  the  model  inputs. 
Although  the  density  jump  within  the  thin  layer  is  indicated,  the  sound  speed 
change  is  not  readily  apparent  at  that  depth.  For  the  highest  frequency 
allowed  in  the  bandlimited  time  responses,  the  thin  layer  is  less  than  one 
wavelength  thick.  Higher  frequencies  are  required  to  obtain  better  resolution 
in  depth. 

The  time  responses  obtained  using  the  same  Kaiser  window  filter  described 
earlier  are  presented  in  figure  17.  Now  the  continuous  reflections  from  the 
thick  layer  are  more  easily  observed.  With  these  time  responses  used  in  the 
inversion  algorithm,  the  reconstructed  results  shown  in  figure  18  are 
obtained.  Except  for  the  poor  resolution  of  n(z)  near  the  thin  upper  layer, 
the  reconstructions  agree  well  with  the  model  inputs,  particularly  for  the 
density  values.  The  value  of  c0/c{z)  for  z  >  300  m  is  slightly  in  error. 
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Figure  13.  (a)  Amplitude  and  (b)  Phase  of  the  Passband  Frequency  Response  of 

the  Reflection  Coefficient  for  Model  3,  Grazing  Angle  e0  =  60’ 
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Figure  17.  Bandlimited  (0-64  Hz)  Time  Responses  of  the  Reflection 
Coefficients  for  Model  3  at  Grazing  Angles  (a)  e0  =  60'  and  (b)  e0  *  90*, 

Kaiser  Window  Filter 


Figure  18.  Reconstructed  Profiles  for  Model  3  Using  the  Time  Responses  in 
Figure  17  Showing  (a)  Normalized  Density  and  (b)  Refractive  Index 
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7.  SUMMARY 


The  inversion  method  of  Candel  et  al.4,5  has  been  implemented  with  a 
view  to  recovering  the  acoustic  properties  of  a  layered  ocean  bottom.  For  the 
scattering  of  acoustic  plane  waves,  the  method  permits  the  reconstruction  of 
both  the  density  and  the  sound  speed  profiles  via  numerical  integration  of  a 
system  of  four  first-order  differential  equations.  Reflection  data  for  two 
distinct  grazing  angles  are  required.  To  test  the  code,  noise-free  reflection 
data  in  the  form  of  bandlimited  impulse  responses  were  generated  synthetically 
for  three  nonabsorbing  inhomogeneous  models.  The  results  of  the  numerical 
inversions  are  in  good  agreement  with  the  original  models.  However,  the 
effects  of  absorption  and  noisy  reflection  data  require  further  investigation. 
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Appendix 


RICATTI  EQUATION  FOR  U/D 


While  the  local  fields  U  and  0  satisfy  the  coupled  system  in  equation 
(20),  the  reflection  coefficient  R  =  U/D  satisfies  a  nonlinear  differential 
equation.  It  is  straightforward  to  show  that 

R'  =  (U/D)'  =  U'/D  -  UD'/d2  =  -2ikzR  +  g(l  -  r2),  (A-l) 

where  U'  and  O'  were  eliminated  via  the  system  in  equation  (20).  Equation 
(A-l)  is  of  the  Ricatti  type.  From  initial  conditions  (22)  and  (23),  it 
follows  that 


R(H)  =  0.  (A-2) 

Recalling  the  jump  conditions  (28)  and  (29)  on  the  U  and  D  fields  at 
discontinuities  in  the  medium,  it  can  be  shown  that  the  ratio  U/D  undergoes  a 
jump  condition  given  by 

R_  -  [d  -  W_)  +  (1  +  VY_)R+][(1  +  Y+/Y_)  +  (1  -  Y+/Y_)R+]-1  .  (A-3) 


This  Ricatti  equation  is  readily  converted  into  the  integral  form, 

z  y  s 

R(z)  =  /  g(s)  [1  -  R(sr]  exp  [2i/k  (t)  dt]  ds,  (A-4) 

H  z  Z 


where  the  limits  of  integration  were  adjusted  in  the  same  way  as  that  which 
preceded  equation  (37). 

The  success  of  the  approximation  given  in  equation  (37)  appears  to  rely 
on  | B |  «  | A  |.  This  suggests  that  equation  (A-4)  can  be  solved  by  iteration 
whenever  U/D  is  small.  In  the  zeroth  approximation,  we  substitute  R(°)  =  0 
into  the  integrand  of  equation  (A-4)  to  obtain  for  the  next  iterate,  when 
evaluated  at  z  =  0,  the  result, 


R(l) 


(0)  =  - 


H  S 

/  g(s)  exp[2i  ;  k  (t)  dt]  ds, 
0  0  z 


(A— 5 ) 


which  is  in  agreement  with  equations  (37)  and  (47). 
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